#######################################################################################
#######################################################################################
####Non-state Atrocities in Capital Cities - Urbanization Levels (2009) Global Map ####
#######################################################################################
#######################################################################################
library(foreign)
library(doBy)
library(plyr)
library(ggplot2)
library(ggthemes)
library(maps)
library(maptools)
library(ggmap)
library(mapdata)
library(RgoogleMaps)

# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")

###Create data for plotting
#Collapse urbanization levels by coordinates
values<- summaryBy(lagurban ~ latitude + longitude, data=main.data, FUN=c("mean"), na.rm=TRUE)
names(values) <- c("latitude", "longitude", "lagurban")
#Create logged values
values$lagurban <- ifelse(values$lagurban>0, (values$lagurban+1), 0)
#Normalize values by relations to to the grid cell with the highest level of urbanization in the sample
values$lagurban <-values$lagurban/max(values$lagurban)
#Create a color spectrum for the map
col.s <- rgb(0, 0, 1, alpha=1*values$lagurban, maxColorValue = 1)
#Create color spectrum for the legend
c0 <- rgb(1, 1, 1, maxColorValue = 1)
c1 <- rgb(0, 0, 1, alpha=0.25, maxColorValue = 1)
c2 <- rgb(0, 0, 1, alpha=0.50, maxColorValue = 1)
c3 <- rgb(0, 0, 1, alpha=0.70, maxColorValue = 1)
colvec <- rbind(c0, c1, c2, c3)
#Create categories for the legend
leg.txt <- c("No Urbanization", "Low Urban Density", "Medium Urban Density", "High Urban Density")
###Plot the map
pdf("urbmap.pdf")
map("worldHires", col="white", fill=TRUE,xlim=c(-180, 180), ylim=c(-58, 90))
points(values$longitude, values$latitude, pch=14, col=col.s, cex=0.2)
legend("bottomleft",title = "Change in Urbanization Density", 
       legend = leg.txt, 
       fill=colvec,
       cex = 0.56,
       bty = "n")
dev.off()
